Sample mixups in Attie eQTL data

There were originally 554 mice. There were 3 mice with failed genotypes. We omitted the genotypes for another 7 mice, as the genotype data were bad. (They showed a high rate of apparent genotyping errors, an unusually large proportion of homozygous gneotypes, and many apparent crossovers.) That leaves 544 mice with genotype data.

The genotype data were cleaned to remove markers not segregating in the cross, to correctly align the BTBR and B6 alleles, and to remove some badly behaved markers (data not shown). Remaining were 2060 markers, including 20 on the X chromosome.

We have gene expression microarray data on six tissues: adipose, gastroc, hypo, islet, kidney, and liver. These were two-color Agilent arrays, with the probes being 60-mers. For each array, one color was a tissue-specific pool and the other was an individual sample. The expression levels were quantified as the “mlratio”, which is approximately the log₁₀ ratio of the individual to the pool, though with some corrections that I don’t fully understand, and with values always between -2 and +2.

There were approximately 500 expression arrays performed for each tissue. A number of outlying arrays were omitted from each group. For hypo, there was an additional batch of 119 poorly behaved arrays that were omitted in later analyses, but they are included here as they contain a number of detectable sample mixups.

There were 493 individuals assayed for adipose expression, after omitting 4 bad arrays. The analysis described above produces a 493 × 527 of median correlations, with the rows corresponding to the adipose samples.

The following figure contains histograms of the self-self correlations and self-non self correlations (or really medians of correlations). The self-self correlations are mostly large, but there are 5 samples with median self-self correlations < 0.5. The two modes in the self-nonself correlations correspond to mixed-sex and same-sex pairs.

For each row, we pull out the maximum correlation (or really the maximum of the median correlations); in the following, we plot the self-self correlations against these maxima. For most samples, the self-self correlation is the maximum; these samples appear to be aligned correctly. The 5 samples with problems are highlighted in green, and we indicate the sample label and then what we infer to be the correct label.

For the 5 problem samples, the self-self correlation is low but in each case there is another sample for which correlation is high. We infer that Mouse3583 and Mouse3584 were swapped, and then there is a rotated trio, Mouse3187 → Mouse3188 → Mouse3200 → Mouse3187.

Let’s look at detailed scatterplots of the expression data for the sample swaps. We pull out the probes that show correlation > 0.6 between adipose and at least one other tissue, and then make a scatter plot of a sample from one tissue against a sample from the second tissue. Points correspond to probes, with the sizes of the points indicating the probes' correlations, across individuals, for that pair of tissues.

First, Mouse3583 and Mouse3584. You can see that, in all tissues, Mouse3583 is not correlated with itself but is correlated with Mouse3584, and vice versa.

Now the rotated trio, Mouse3187 → Mouse3188 → Mouse3200 → Mouse3187. Note that Mouse3188 was swapped with Mouse3179 in hypo (see below). Mouse3187 in hypo seems correct but poorly behaved. Mouse3187 and Mouse3188 were not assayed for kidney expression.

There were 496 individuals assayed for gastroc expression, after omitting 2 bad arrays.

We jump to the plot of the self-self correlations versus the maximum correlation. Again, for most samples, the self-self correlation is the maximum; these samples appear to be aligned correctly. The 2 samples with problems are highlighted in green, and we indicate the sample label and then what we infer to be the correct label.

There is an additional potential problem highlighted in orange, having low self-self correlation, and being slightly more correlated to another sample than to itself. But this is Mouse3188, which was involved in sample mix-ups in both adipose and hypo. We believe it’s correctly labeled in gastroc.

The 2 problem samples are inferred to be a sample swap, Mouse3655 ↔ Mouse3659.

Let’s look at detailed scatterplots of the expression data for the sample swap, Mouse3655 and Mouse3659. You can see that, in all tissues, Mouse3655 is not correlated with itself but is correlated with Mouse3659, and vice versa.

There were 493 individuals assayed for hypo expression, after omitting 1 bad array.

Here is the plot of the self-self correlations versus the maximum correlation. Again, for most samples, the self-self correlation is the maximum; these samples appear to be aligned correctly. The 17 samples with problems are highlighted in green, but there are too many to attach labels.

Here are the summary statistics for the problems.

Mouse3589 and Mouse3590 have reasonably high self-self correlation (0.7) but are each much more highly correlated to the other (0.9).

We infer a set of nine sample swaps:

There is one wrinkle here: Mouse3381 was not assayed for hypo expression, but really it was Mouse3382 that was not assayed; Mouse3381 was assayed but was labeled Mouse3382.

Let’s look at detailed scatterplots of the expression data for all of these sample swaps.

First, Mouse3179 ↔ Mouse3188. Recall that Mouse3188 was part of a rotated trio in adipose (Mouse3187 → Mouse3188 → Mouse3200 → Mouse3187). These are not too pretty, but it is a clear swap Mouse3179 ↔ Mouse3188 in hypo.

Mouse3208 ↔ Mouse3210: These are rather messy, but it is still a clear sample swap.

Mouse3347 ↔ Mouse3348: Again, these are rather messy, but still a clear sample swap.

Mouse3367 ↔ Mouse3369

Mouse3381 ↔ Mouse3382 (but note that Mouse3381 was not intended to be done in hypo).

Mouse3449 ↔ Mouse3451

Mouse3452 ↔ Mouse3454

Mouse3589 ↔ Mouse3590

Mouse3592 ↔ Mouse3594

There were 498 individuals assayed for islet expression, after omitting 1 bad array.

Here is the plot of the self-self correlations versus the maximum correlation. Again, for most samples, the self-self correlation is the maximum; these samples appear to be aligned correctly. The 3 samples with problems are highlighted in green, and we indicate the sample label and then what we infer to be the correct label.

There is an additional potential problem highlighted in orange, having low self-self correlation, and being slightly more correlated to another sample than to itself, but this is Mouse3188 again, which was involved in sample mix-ups in both adipose and hypo. We think it is correctly labeled in islet.

The 3 problem samples are inferred to be a sample swap, Mouse3598 ↔ Mouse3599, plus a sample duplicate: Mouse3295 was run correctly but also with the label Mouse3296 (more on this below).

Let’s look at detailed scatterplots of the expression data for the sample mixups.

First, the sample swap, Mouse3598 ↔ Mouse3599. You can see that, in all tissues, Mouse3598 is not correlated with itself but is correlated with Mouse3599, and vice versa.

Now we turn to the apparent duplicate, the Mouse3296 sample really being Mouse3295. Mouse3296 in islet is not correlated with itself but is correlated with Mouse3295. Mouse3295 in islet looks fine: It’s correlated with itself and not correlated Mouse3296 in the other tissues. We’ll revisit this below in the section on Sample duplicates.

There were 481 individuals assayed for kidney expression, after omitting 1 bad array. Note that there were 27 samples that were assayed for kidney expression but no other tissue.

Here is the plot of the self-self correlations versus the maximum correlation. Again, for most samples, the self-self correlation is the maximum; these samples appear to be aligned correctly. The 2 samples with problems are highlighted in green, and we indicate the sample label and then what we infer to be the correct label.

There is an additional potential problem highlighted in orange (Mouse3484), being a bit more correlated to another sample (Mouse3503) than to itself. These two arrays are oddly behaved and will be investigated further below.

The 2 problem samples are inferred to be a sample swap, Mouse3510 ↔ Mouse3523.

Let’s look at detailed scatterplots of the expression data for the sample mixups.

First, the sample swap, Mouse3510 ↔ Mouse3523. You can see that, in all tissues, Mouse3510 is not correlated with itself but is correlated with Mouse3523, and vice versa.

Now we show scatterplots for the odd pair, Mouse3484 and 3503. They all seem correlated with each other. We’ll revisit this below in the section on Sample duplicates.

There were 490 individuals assayed for liver expression, after omitting 1 bad array.

Here is the plot of the self-self correlations versus the maximum correlation. The 2 samples with problems are highlighted in green, and we indicate the sample label and then what we infer to be the correct label.

The 2 problem samples are inferred to be a sample swap, Mouse3142 ↔ Mouse3143 (though Mouse3143 was not intended to be assayed for liver), and a sample duplicate: Mouse3136 was assayed correctly and also with the label Mouse3141 (more below).

Let’s look at detailed scatterplots of the expression data for the sample mixups.

First, the sample swap, Mouse3142 ↔ Mouse3143, but with Mouse3143 not intended to be assayed in liver. You can see that, in all tissues, Mouse3142 is not correlated with itself but is correlated with Mouse3143.

Now we turn to the apparent duplicate, the Mouse3141 sample really being Mouse3136. We’ll revisit this in the next section, on Sample duplicates.

Above, we had identified several apparent sample duplicates (which might be called “unintended technical replicates”). We now look at these in more detail.

Consideration of all probes on an array will not be informative, as most genes are not expressed in a tissue and so the measurements are just noise. To identify probes that with signal, we pull out those probes with correlation > 0.6 against at least one other tissue.

For each tissue, there are non-duplicate pairs with correlations that reach 0.9, while each of islet, kidney and liver have a pair of samples with correlation > 0.98.

Here are the top 5 between-tissue correlations for each tissue.

In islet, Mouse3296 is a duplicate of Mouse3295.

In liver, Mouse3141 is a duplicate of Mouse3136.

In kidney, Mouse3484 and Mouse3503 remain a bit of a conundrum, since they each are correlated with the other in all tissues.

This feature is specific to the two kidney samples. For example, consider the between-individual correlations for all 15 tissue pairs. Here, we’ll grab the subset of 76 probes with correlation > 0.6 for all tissue pairs.

Here are plots of the 15 between-tissue correlations for these two mice. The pairs involving kidney are in red.

I still don’t know what to make of this. It’s interesting to note that the red points are at about the same height in all four figures, while the blue points are high for 3484:3484 and 3503:3503 but low for 3484:3503 and 3503:3484.

Why should this pair of kidney samples have such similar expression? I’m inclined to think that a mixture of the two RNAs were placed on both arrays.

I will omit these two arrays from further analyses.

Here’s a summary of the arrays omitted and found do be in mislabeled, for each tissue.

I now omit the relevant arrays, combine the duplicates in liver and islet (by simply taking the average of the mlratios for each probe) and then re-name the arrays as inferred.

Note that the points are colored based on the results for the six tissues combined. For each tissue there are numerous green points (with the combined data, viewed as “fixable”) that have high minimum mismatches. These are cases where the inferred sample was not run in that particular tissue but was run in others: For the most part, these are samples that were run only in kidney or run in other tissues but not kidney. Hypo is a bit messy, though, with less-clear separation between “fixable” and “not fixable”.

Here is a table of these discrepancies. The numbers are the minimal percent mismatches between observed and inferred eQTL genotypes, by tissue and then with all tissues combined. The “Note” column explains why some tissues show a discrepancy; most are due to missing expression data.

There are two additional samples that deserve to be highlighted, which are two blue points (inferred to be correctly labeled) that are off the diagonal line in one tissue:

In both of these cases, I’m inclined to conclude that the expression array data may be of low quality, but the samples are correctly labeled.

The re-alignment of the DNA samples requires that we go back to the inferred sample duplicates and relabel some of them. In particular, the sample for Mouse3269 was found to be a duplicate of the sample for Mouse3259, but the sample for Mouse3259 was actually for Mouse3332, and so both are really Mouse3332. Similarly, the sample for Mouse3354 was a duplicate of the sample for Mouse3353, but the sample for Mouse3353 was actually Mouse3352, and so both are really Mouse3352.

Here’s a summary of the problems observed in the genotype data.

I originally came to these observations by identifying a set of samples for which sex did not match the X chromosome genotype. It would be good to revisit that issue: with the genotypes realigned, do the sexes now match what was observed in the X chromosome data?

We first grab the set of DNAs for which we are now sure of the appropriate labels. (There are 519 such.)

We want to compare the inferred sex from the X chromosome genotypes for the samples to the documented sex (attached to the correct label).

Of the 519 DNA samples that are correct or can be corrected, 23 showed a mismatch between sex and X chromosome genotypes. With the corrected sample labels, there are 3 discrepancies:

But all three of these mice are homozygous (or hemizygous) BTBR across the entire X chromosome, and so the X chromosome genotypes are compatible with either sex. Thus, there are no real discrepancies between the X chromosome genotypes and the sex of the mice, once the labels for the DNA samples have been corrected.

We can also look at the expression of the Xist gene (which should be high in females and 0 in males) and of Y chromosome genes (which should be >0 in males and 0 in females).

In the following, I plot the average expression of 3 Y chromosome genes against the expression of Xist, with points colored by the sex of the sample (blue is male; red is female). (The expression arrays contained 21 probes for genes on the Y chromosome, but only a portion of them showed clear sex differences in expression in all tissues, and so I focused on these.)

Hypo is a bit of a mess, with a bunch of females having unusually low Xist expression and a bunch of males having unusually low Y chromosome expression; this is almost entirely explained by a bad batch of arrays that we ultimately removed.

The other five tissues look great. There is one case of a male with relatively high Xist expression in adipose, but the sex assignments of the expression arrays, after the fixing labels, appears to be correct.

The following is the plot for hypo, without the bad batch of arrays:

First, I’ll look at insulin levels at 10 weeks; we actually have three different measures of insulin at 10 weeks, and we also have insulin at 4, 6 and 8 weeks, but I’ll just look at one of the 10 week measures. I’ll use the log scale and will include sex as an interactive covariate (that is, allowing QTL to have different effect in the two sexes).

In the original data, there were 539 mice; in the corrected data, we are left with 519 mice. The following are the LOD curves, with those for the original data in red and those for the corrected data in blue.

With the original data, there were 3 chromosomes having a LOD score ≥ 5; with the corrected data, there are 7 chromosomes with LOD ≥ 5.

Here is a summary table of the peak LOD on chromosomes with inferred QTL.

The LOD scores on chromosomes 2, 5, 6, 7, and 9 went up substantially. The peak locations are largely unchanged.

For most mice, it was recorded whether they had an agouti coat or not. This is due to a single gene (a, non-agouti) on chromosome 2 at 154.6-154.8 Mbp.

In both the original and corrected data, agouti coat color maps strongly to chromosome 2. The LOD score increased from 64 to 110 with the correction of the sample mixups. In both cases, the peak is at 154.4 Mbp, right next to the true gene.

If we grab the mice with good genotype data at the agouti locus, we obtain the following relationship between agouti genotype and phenotype for the original data (with alleles B = B6 and R = BTBR):

With the corrected data, we get the following table.

With the original data, there are 47 mismatches out of 498 mice. With the corrected data, there are 7 mismatches out of 480 mice.

So there are still 7 mismatches remaining, but the results are considerably improved.

Here are the actual mice with mismatching agouti genotype/phenotype:

The genotypes for the mismatching individuals seem clear and are not likely genotyping errors. Here are the genotypes on chromosome 2 for these seven mice, with the agouti locus indicated by a green vertical line. The top five have tan coats and should be BR (gray points) or RR (black points) at the agouti locus, but are all BB (white points). The bottom two have black coats and should be BB (white points) at the agouti locus but are BR (gray points).

For most mice, we also have information on whether they had a tufted coat or not. This is due to a single gene on chromosome 17, though it seems that the actual gene is not yet known.

In both the original and corrected data, tufted coat maps strongly to chromosome 17, though the LOD score increased from 64 to 107 with the correction of the sample mixups. In both cases, the peak is at marker rs3700924, at 27.2 Mbp.

If we grab the mice with good genotype data at the tufted locus, we obtain the following relationship between tufted genotype and phenotype for the original data (with alleles B = B6 and R = BTBR):

With the corrected data, we get the following table.

With the original data, there are 37 mismatches out of 538 mice. With the corrected data, there are 4 mismatches out of 519 mice.

So there are still 4 mismatches remaining, but, as with agouti coat, the results are considerably improved after the correction for the sample mixups.

Here are the actual mice with mismatching tufted genotype/phenotype:

The genotypes for the mismatching individuals seem clear and are not likely genotyping errors. Here are the genotypes on chromosome 17 for these four mice, with the tufted locus indicated by a green vertical line. These all have non-tufted coats and should be BB (white gray) or BR (gray points) but are all RR (black points) at the tufted locus.

Finally, we will map the expression of genes in all six tissues, with both the original data and the corrected data. We will focus on probes for which we were able to infer a genomic location, and on an autosome or the X chromosome (i.e., not the Y or mitochondria). There are 37797 such probes.

For all probes, we transform the expression measures to normal quantiles, and we include sex as an interactive covariate in the QTL analyses.

For the hypo expression data, we omitted a batch of 119 arrays that were found to be bad. For each transcript in each tissue, we used Haley-Knott regression to calculate LOD curves across the chromosome. For each transcript, we find the maximum LOD score on each chromosome and count the number of chromosomes that have LOD ≥ 5. A transcript is said to have a local-eQTL if the LOD score on its chromosome is ≥ 5 and the genomic location of the transcript is within the 2-LOD support interval for that eQTL; otherwise, a peak above 5 on that chromosome is considered a trans-eQTL.

Here is a figure with the numbers of local- and trans-eQTL for each tissue, with the original data and with the corrected data.

Across all tissues, there were 16946 local-eQTL and 69332 trans-eQTL with the original data. After correcting the sample mixups, these numbers increased to 18183 and 95152 for local-eQTL and trans-eQTL, respectively. These are increases of 7% and 37%, respectively.

In response to a question from a reviewer, let’s repeat this, looking for eQTL with LOD ≥ 10.

Across all tissues, there were 9527 local-eQTL and 12262 trans-eQTL with the original data. After correcting the sample mixups, these numbers increased to 10339 and 17470 for local-eQTL and trans-eQTL, respectively. These are increases of 9% and 42%, respectively.